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SUMMARY 

This method is investigated for semi-infinite multiple -slab config- 
urations of arbitrary width, composition, and source distribution. 

Isotropic scattering in the laboratory system Is assumed. 

Isotropic scattering implies that the fraction of neutrons scattered 
in the i "* 3 * 1 volume element or subregion that will make their next colli- 
sion in the volume element or subregion is the same for all collisions. 
These so-called "transfer probabilities" between subregions are calcu- 
lated and used to obtain successive-collision densities from which the 
flux an d transmission probabilities directly follow. 

For a thick slab with little or no absorption, a successive- 
collisions technique proves impractical because an unreasonably large 
number of collisions must be followed in order to obtain the flux. Here 
the appropriate integral equation is converted into a set of linear 
simultaneous algebraic equations that are solved for the average total 
flux in each subregion. 

When ordinary diffusion theory applies with satisfactory precision 
in a portion of the multiple-slab configuration, the problem is solved 
by ordinary diffusion theory, but the flux is plotted only in the region 
of validity. The angular distribution of neutrons entering the remaining 
portion is determined from the known diffusion flux and the remaining 
region is solved by higher order theory. 

Several procedures for applying the numerical method are presented 
and discussed. To illustrate the calculational procedure, a symmetrical 
slab in vacuum is worked by the numerical, Monte Carlo, and P 3 spheri- 
cal harmonics methods. In addition, an unsymmetrical double-slab problem 
is solved by the numerical and Monte Carlo methods. The numerical 
approach proved faster and more accurate in these examples. Adaptation 
of the method to anisotropic scattering in slabs is indicated, although 
no example is included in this paper. 
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INTRODUCTION 

In order to solve a multiple -slab configuration of arbitrary width, 
composition, and source distribution for the monoenergetic neutron flux 
distribution and the number of transmissions, recourse to a number of 
powerful methods is available. Wherever diffusion theory does not apply, 
the spherical harmonics (ref. l), Monte Carlo (ref. 2), and S (ref. 3) 
methods can be used. 1 

The present discussion concerns a numerical method that accurately 
approximates the exact solution to many slab problems with a reasonable 
amount of calculation. Isotropic scattering of monoenergetic neutrons 
in the laboratory system is assumed. The numerical method is essentially 
a procedure for solving the integral equation for the total flux given 
in reference 1 (p. 78). 

The numerical method was worked out in principle by DeMarcus and 
Nelson (ref. 4) and was obtained by the author as a natural consequence 
of working slab problems by Monte Carlo. This method is relatively easy 
to apply when at least some absorption is present and/or for thin slabs 
where an appreciable leakage occurs - in other words, where too many 
collisions need not be followed in order to obtain the flux; these are 
the conditions under which diffusion theory is generally inadequate. 

For the other extreme, a thick slab with little or no absorption, a set 
of simultaneous algebraic equations can be solved for the average total 
flux in each subregion. When the diffusion- theory solution is known to 
apply with satisfactory precision in a portion of the slab array, this 
information can be utilized to great advantage in simplifying the numeri- 
cal calculation. Also, it proves feasible in a number of cases to treat 
a multiple-slab configuration as essentially a single-slab problem with 
the resulting simplification. 

For purposes of illustration, a simple slab containing a centrally 
located isotropic source plane and surrounded by vacuum was chosen. A 
scattering probability of 8/ 10 per collision and a slab half-width of 
0. 781 mean free paths were also chosen. 

The integrals for the probabilities of transmission with no colli- 
sion and one collision were analytically formulated and evaluated. The 
numerical -method flux and the number of neutrons transmitted with 
exactly K collisions, K = 0,1,2, . . . 19, were obtained and compared 
with a Monte Carlo calculation. The numerical -method flux was also 
compared with the P^ spherical harmonics aid P 1 diffusion-theory 
fluxes. 

Values of the slab half -width were successively changed to l/s and 
5 total mean free paths, and the numerical -method flux was obtained and 
compared with diffusion theory for these cases. 
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As a second example, a configuration in vacuum consisting of two 
contiguous slats of different absorption and scattering properties was 
solved. An isotropic source plane was assumed at the left extreme 
boundary. Distributed volume sources were avoided to facilitate the 
Monte Carlo calculation. Twenty-eight collisions were followed in the 
Monte Carlo and numerical successive-collisions treatment. The total 
flux and the number of transmissions after exactly K collisions were 
recorded and compared. The total flux was also obtained numerically by 
solving the appropriate set of simultaneous algebraic equations discussed 
in the text. 


SYMBOLS 

regions or slabs A,B,C, . * . 
arbitrary constants 

widths in centimeters of slabs I, II, III, . . . 
or slabs A,B,C, . . . 

a 0, a l, . . . bQ,b^ nuclear constants of region A and B (defined in 

appendix F) 

F,F(x,|i) transport flux 

Fj(x) coefficient of the Legendre polynomial of \i 

in series expansion of transport flux 

matrices defined by equations (3) and (12) 

number of spherical zones of width dp represent- 
ing equal areas on a hemisphere of unit radius 

XX y XX -j ^ 

matrices defined by equations (3), (12), and (13) 

H 0,2 ;H 1,2 

I q-by-q- unit matrix 

i=l,2,3, . . . index denoting specific subregion (subslab) of 

given slab (subregions are numbered consecutively 
from left to right) 

j=l, 2, 3, . . . index denoting specific subregion (subslab) of 

given slab 

K number of collisions a neutron has suffered 


G ’\,0’%,q 2 

g 


A, B, 0, • • • 
A v ,B v ;v=l,2,3,4 
a,b, c, . . . 
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£L 

L' 


L 


r 


l,V 


N S ,I^ 

n|(i) 


n(x 0 ,n)dM. 


nS (i) 


n 


T 


T 

"K 


subregion width in total mean free paths 

distance in mean free paths from point lying any- 
where in subregion i to center of subregion j; 
neutron is assumed to scatter in i and then 
collide in j 

distance in mean free paths from center of subregion 
i to center of subregion j , where 
r = |j-i | = 1,2 . . . 

distances in mean free path units defined in connec- 
tion with equation (A3) and sketch (c), page 26 

distance in mean free paths from boundary of one 
subregion to center of another or conversely 

matrices defined by equations (5) and (3) 

number of neutrons scattered per second from 
collision in subregior ij in particular, n^(i) 
represents primary volume source in i 

number of neutrons per square centimeter per second 
arriving at Xq and laving direction cosines in 
dji about fj. 

total number of neutrons per second scattered in i 

for all collisions . v 

total number of neutrons transmitted after one or 
more collisions 

number of neutrons tranEmitted in exactly K 
collisions 


P 

P(AL p £L I:[ ) 
P(x 1 , x)dx 


matrix defined by equation (15) 

transfer probability between unequal subregions as 
measured in units of mean free path 

kernel of equation (10), given by equation (A2) 


p i,j ° r 

PfXi.XjjAXj 


probability that a neutron will make its next colli- 
sion in subregion j after being isotropically 
scattered in subregior i 


P l(n) 



Legendre polynomial 


* 
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P r or 
P(L r ,AL) 


pLR 


P T (Z) or 
P T (x’,x) 



K 


p(q’ ,q)dq 


P r ; r=l,2, . . . 

q. 

R 

r 

S 

S(L' )dL' 

AX 

X i’ X j 
x 0> X L* X R 

P V >Y V > VJ 
\ ,X 1' X 2 [ 


transfer probability between subregions of equal 
width, as measured in units of total mean free 
path ^eq. (A8a) 

probability that next collision will occur in sub- 
region in which neutron scattered 

probability that a neutron will make its first 
collision to left of source plane and be trans- 
mitted through right-hand boundary without a 
second collision 

probability that a neutron scattered or born iso- 
tropically at x' will be transmitted through a 
boundary Z mean free paths away with no inter- 
vening collisions 

probability of a source neutron being transmitted 
through extreme left or right boundaries, respec- 
tively, in exactly K collisions 

given a neutron having direction cosine q' and 
undergoing scattering collision that may or may 
not be isotropic; this expression represents 
probability of finding neutron velocity in direc- 
tion dq about q after collision 

transfer probability from source plane at boundary 
of a subregion to within r th adjacent subregion 

total number of subregions in given configuration 

random number 

integer 1,2,3, . . . equaling j-i 
source 

source in infinitesimal strip of width dl' 
subregion width in centimeters 

x coordinates of centers of i th and j th subregions 

x coordinates of interface and extreme left and 
right boundaries, respectively 


defined by equation (F3) 
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&(x) 


Dirac delta function 


^,1 


probability that a neutron with direction cosine 
will pass through a j. ingle subregion without 
collision 


Kronecker delta = 


0,h ? 
l ; h ^ 


1 

1 


^h 


€ error per collision resulting from assumption that 

all scattering collisions within a subregion take 
place at center 

U(V\) total mean free path anc. scattering probability per 

collision 


P 

I 

P 

Z 

0 

<p(i) 

Q. 

aa 


equals cosine of £ 
angle between ft and the x axis 
radial distance to next collision point 
equals l/A 

1-by-q matrix for flux in each subregion 
average total flux in sibregion i 
unit vector in direction of neutron velocity 
element of solid angle about Q 


Subscripts and superscripts: 


h 

L 

R 

s 

T 

I, II 


specific solid angle in direction space 

left extreme boundary 

right extreme boundary 

scattering 

transmitted 

specific slab, medium, cr region 


t?d 




AmLYSIS 

The assumption of isotropic scattering in the laboratory system is 
a good approximation for neutron collisions \sith heavy nuclei. Isotropic 
scattering implies that the direction of neutron travel after collision 
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is unrelated to the direction of travel before collision. Thus, the 
probability P., * that neutrons scattering in volume element i will 

make their next collision in volume element j is the same for all 
collisions. For a one -dimensional or semi-infinite slab, i and j 
refer to subslabs or subregions that hereafter will be taken to be equal 
in width in a given slab. This implies that Pq^j = Pj,i* Furthermore, 

if these subregions are numbered consecutively, P^j = p i+l,j+l* This 
means, for example, that P-^ 3 = ^ 3,5 = ^6,4 = ^8,6' ^ so fort ' h * ' I ' he 
transfer probability, henceforth called P r , therefore depends on a 
single index r and is characteristic of slab geometry. Thus, given 
that a neutron is scattere in any subregion, Pq is the probability 
that it will make its next collision in the same subregion, Pq in 
the adjacent subregion, Pg in the second adjacent subregion, and so 
forth. The P r 's do not vary from collision to collision and can there- 
fore be calculated once and used to obtain the (K+l) s ^ collision density 
distribution from the K th . For a slab with q equal subregions, there 
are q different transfer probabilities. 

For problems involving spherical symmetry, the sphere would be 
divided into q spherical shells of equal width. Here, however, the 
transfer probability P^ j from shell i to shell j is not equal to 
that from j to i, P^ j. Also, / ^j,j* '® 1US > there are q 

different transfer probabilities as opposed to the q probabilities of 
the slab case. 

In cylindrical geometry there would also be q^ Pq^j' s between 
concentric cylindrical shells. This paper considers only the slab case. 


Single -Slab Analysis 


Consider a typical slab of width "a" centimeters divided into q 
subregions of equal width, numbered from left to right by the index I. 


A s be the total and scattering mean free paths of the 

AT, total mean free 
to be the number of neutrons scat- 


Let A and 

given medium and let the width of each subregion be 
paths. Furthermore, define n|(i) 
tered per second from the K th collision in subregion i. Then n 
for example, represents the given initial volume source in i, and 


(i)> 


iUl j JL ^ o— 

n s ( i) represents the number of first-collision scatterings in 
AT. has been chosen sufficiently small. 


1. 


If 


&L = 


AK 


T = qX 
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then the sources within each subregion can te lumped at its center with 
negligible error (appendix A). The transfei probability is then computed 
from the center of the i th subregion to witi in the i+r th subregion, after 
which the transferred neutrons are in turn Jumped at the center of the 
latter subregion. 

The transfer probability P r is a function only of the distance 
L r between the centers of the two subregions involved and of the sub- 

region width AL, all in units of total mean free path. Although 
L r = rAL, r = 1,2, . . (q-l), P p is denoted by P(L r ,AL) as formu- 

lated and evaluated in appendix A. 

Initial and subsequent collision densit ies arising from distributed 
isotropic volume sources . - The specified volume sources are assumed to 
be lumped into isotropic source planes located at the centers of the sub- 
regions. After , ^q-l are calculated from equations (AQa) 

and (A9), the number of (K+l) st collision scattering hits in subregion i 
may be obtained by summing the contributions to i from the previous or 
K th collision in all of the subregions. Mat nematic ally stated, 



where the first summation represents the input into i from subregions 
to the left of i and the second summation represents the input from 
subregions to the right of i, including i itself. The quantity 
A/A s gives the fraction of the (K+l) st collision hits that scatter. «" 

Setting K = 0 in equation (l) gives the first-collision distribution 
due to the distributed primary volume sources n^(i). Successive- 

collision distributions are obtained by successively increasing K by 
1 unit and evaluating n^ +1 (i) from equatioi (l). The total number of 

scattering hits for all collisions is given by 


n S (i) = 2 n|(i) 

K=1 


( 2 ) 


A sufficient number of collisions are followed until n s (i) does not 
change appreciably. 


First-collision distribution due to isosropic source plane at inter- 
face between subregions . - If the source plane is located at an Inter- 
face between subregions, another set of transfer probabilities p r can 
be calculated from equation (A8) as before, out with 


t 




can 
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L ' = (r - -^AL r = 1,2, . . . (A8b) 

instead of 

L' = r AL r = 1,2, . . . (A8a) 

Thus, represents the fraction of the source-plane neutrons that will 

collide in the first adjacent subregion, p^ in the second adjacent sub- 
region, p 3 in the third, and so forth. Once the initial-collision dis- 
tribution has been determined, successive-collision distributions follow 
from equation (l). 

If a reasonable amount of absorption is present and/or the slab is 
sufficiently thin so that a large number of collisions are not required 
in order for the number of scattering hits to converge, then equation 
(2) may be utilized for machine calculation. For cases where conver- 
gence is slow, a collision-by-collision technique is impractical, so a 
matrix treatment, suggested in reference 4, is used. 
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Equation (l) can then he written in matrix form as 

l£ +1 = h£h K = 0,1,2, ... (4) 

The symbol Nq defines the 1-by-q matrix for the given original volume 
sources. The matrix to be found is N s , a 1-by-q matrix whose elements 
n s (i) represent the total number of scattering hits in subregion i for 
all collisions (eq. (2)). 

Matrix N s = [n s (l), n s (2), . . ., n s (q)] 

From equation (4) it follows that 

*Li - n o hK+1 


or 


oo 



where I is a q-by-q - unit matrix. Therefore, 

N S ( I-H) = NqH (5) 

The quantities H, I-H, and Nq are known matrices given by equation 
(3). Hence, equation (5) represents q linear simultaneous algebraic 
equations for the unknown n s (i). If n s (i) has been found, the average 
total flux in each subregion is computed as 

cp(i) = n s (i) g i = 1,2, . . ., q (6) 

where AX is the width of subregion i in centimeters. In some cases 
it will be more economical to utilize equation (5) to obtain the con- 
verged flux, while in others a collision-by- collision technique as 
implied by equation (4) will be preferable. 


Integral equation for total flux . - It vill now be shown that equa- 
tion (5) corresponds to numerically solving an integral equation for the 
total flux. 


Denote the x coordinate of the center of the ith subregion in 

which the neutron scatters by x'- . Let X; be the x coordinate of the 
th « 

center of the j subregion to which the neucrons are transferred for 
their next collision. Replace the transfer probabilities in equation (l) 


* 


(V T- 



CY-2 back 


11 


"by their equivalent P(X|,Xj)AXj and replace n^ +1 (i) ^y n^ + ^(Xj)Z*Xj. 
Summing hoth sides of equation (l) over-all collisions gives 


Z - 

K=0 



nJ(Xp^ P(xpX J )AX j 


+ 


X, X 4.(n)^± p(xi,xj)AXj 

i=l K=1 


( 7 ) 


Equation (7) states that the total scattering rate for all colli- 
sions in AXj is the number of first-collision scatterings plus the 
number of subsequent scatterings as represented by the last term* From 
equations (2) and (6), there follows, after dividing through by 


®(Xj) - A V ng(Xi)tfq P(Xi,Xj) + ^- X <5(X!)AX^ P(Xi,Xj) (8) 
1=1 S 1 = 1 


With the subscript i and/or j affixed to hq , P, and 0, equation 
(8) becomes 


Oj = A 


q 

Z 

1=1 






A = i,2. 


q O) 


Equation (9) represents q linear simultaneous algebraic equations. In 
matrix form, it is identically equation (5) with N s replaced by 

s 

Allowing £X.[ -*-0 in equation (8) gives the Fredholm -type integral 
equation for the total flux (ref. l) : 


$(x) = A f nS(x') P(x' ,x)dx' + — f $(x') P(x',x)dx' 

A* A s A* 


(10) 


The kernel P(x',x) is given by equation (A2). The successive-collisions 
technique of equation (4) therefore corresponds to solving equation (10). 
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Multiple-Slab Analysis 


Isotropic 



Sketch (a). - Multiple-slab CDnfiguration. 

Angular distribution of neutrons . - Sketch (a) displays two contig- 
uous slabs. Additional slabs need not be considered because all the 
numerical techniques required are illustrated by this double-slab problem. 
An isotropic source plane is assumed at the interface xq. 

The neutron flux is obtained in the first slab by following neutrons 
from collision to collision within this slab by means of equation (4) or 
by solving the q^ simultaneous algebraic equations represented by 
equation (5). The second slab is entered by means of an angular distri- 
bution n(xQ,|i)dq that is determined from the known flux. The angular 
distribution is derived as equation (B2) and represents the number of 
neutrons per square centimeter per second arriving at xq and having 
direction cosines lying in dp about p. 

The isotropic source plane and n(xQ,p)rlp determine the first- 
collision distribution in the second slab, deutrons are followed therein 
until its flux converges. A new angular distribution is computed for 
reentry into the first slab. These neutrons are followed, giving an 
augmented flux that represents the combined effect of the original neu- 
trons plus those backscattered once from the second slab. The latter 
slab is entered a second time, and so forth. Thus, neutrons are followed 
in one slab at a time, with the other being entered by means of the 
angular distribution. This procedure is not feasible when an excessive 
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number of interface crossings is required to achieve a converged flux 
everywhere. This difficulty does not arise if either slab is a reason- 
able absorber. 

The angular distribution represents an anisotropic source plane at 
the interface. The first-collision distribution of these neutrons is 
computed by dividing the angular hemispherical region - to the right of 
the interface, for example - into a grid of g solid angles. The proper 
number of neutrons is sent into each solid angle by means of n(xQ,p)dp. 
The exponential attenuation law gives the fraction colliding in each 
subregion of the second slab for direction cosine p. Summing over all 
solid angles yields the f irst~collision distribution as expressed by 
equation (B5). Once this is known, subsequent collision distributions 
follow from equation (4), or the final distribution follows from equa- 
tion (5) with in the right-hand term. 

Numerical treatment of anisotropic scattering in slabs . - Aniso- 
tropic scattering in slabs can be treated in an analogous manner by: 

(1) Assuming the anisotropic volume sources for the next collision 
to be lumped into a source plane at the center of a subregion. 

(2) Utilizing an expression p(p ! ,p)dp which gives the probability 
that a neutron having direction cosine p 1 will be found in dp about 
p after undergoing a scattering collision. 

(3) Applying the technique of equations (B5) and (B6) to each source 
plane to determine the number of neutrons that will make their next 
collision in each of the subregions for a given p. 

(4) Repeating steps (l), (2), and (3) for the other q's and pro- 
ceeding to the next collision. The procedure is an order of magnitude 
more involved here because the angular distribution must be recorded at 
the center of each subregion in following successive collisions, whereas 
it was -unnecessary to do this for the isotropic case. 

Reduction of multiple -slab to single-slab configuration . - When the 
number of interface crossings required becomes excessive, the angular - 
distribution approach becomes impractical. Noting that the transfer 
probabilities (eq. (A8)) in a given slab are functions only of the sub- 
region width AL, one attempts to choose ALj = ALjj. The transfer 
probabilities of both slabs become Identical, and a single set of trans- 
fer probabilities applies throughout. A two-slab problem has then been 
reduced to a single-slab problem with the resulting simplification: 

The relation AL-j- = AL-j-j implies that 

a _ b 
q 1 A I q 2 >' I:C 



14 


or 

q 2 = a)* 1 ! = 1 ’ Z >' 5 ’ * • • t 11 ) 

where a and h are the respective widths ii centimeters of slabs I 
and II, A the total mean free path, and q-^ the number of subregions 
in slab I. 

Equation (ll) can be rigorously satisfied only when the given numer- 
ical quantity in parentheses is an integer. Otherwise, a set of values 
should be sought for q n and q^ satisfying (ll) such that: 

(1) Rounding off to the nearest integers alters the given physical 
dimensions of the configuration a negligible amount, 

(2) The values of q n and q^ are not excessively large for cal- 
culation purposes, 

(3) The values of q^ and q^ are sufficiently large that the sub- 
regions have negligible lumping error. 

When these three criteria can be met, it is perhaps best to choose 
the number of subregions in each slab accord! ig to equation (ll). Thus, 
a multiple-slab problem can be worked as a single slab. However, the 
colliding neutrons must be weighted by the scattering probability A/A s 
of the medium in which they collide. Equation (l) must then be modified 
to read: 

"K + l (1 > - (£) £ + S "|( a ) p a-i| 

Lcfcl ct=i J 

= (j£) + Z n K^ a ^ P a- ij 

x b ' [o^l a^i j 

Define matrix ,q identical to matrix G of equation (3 

in the first q^ columns with zero elements In the remaining columns. 
Also define Gq to be identical to G in the last q^ columns with 

zero elements in the others. Defining matrices q and Hq ^ as 


i— 1,2, . . ., q^ 

i = a 1 + 1 > < ! 1 +2, . q 
q. = + q 2 
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H- 


= (XX 


E (' 


H r 


A S/ | 

(X\ lx 


?'s j ^2 


>,2 " (; 

equation (4) becomes 

4 + l = <(=1,0 + H 0, 2 ) S <=1,2 

< 


( 12 ) 


K = 0,1,2, 


(13) 


= [^(l), n®(2), • . n|(q)j 

Following the same procedure used in deriving equation (5) gives 

(14) 


N ( I H l,2^ N 0 H 1,2 


Equation (13) is utilized in a collision-by-collision procedure, 
whereas equation (14) gives n s (i) directly when the q simultaneous 
algebraic equations represented by (14) are solved. 


Transfer probabilities between unequal subregions lying in differ- 
ent media . - When it is impractical to use the angular distribution to 
pass from one slab into the other or to make the AL ! s equal, transfer 
probabilities between unequal subregions lying in different media can be 
used. These transfer probabilities eliminate the need for the angular 
distribution but require many more computer storage locations. They are 
calculated from equation (All). 


Referring to sketch (a), 


P (i) 
U j 


is defined as the transfer probabil- 


ity between subregions i and j lying exclusively in medium I and 
p(2) as the transfer probability for i and j lying exclusively in 


medium II. If i lies in medium I and j in medium II or conversely, 

then bhe transfer probability is written as j and 

respectively. 


A matrix P is defined as 


[(f) 1 p? 1 ! 

V's) 

1 M 11 p - • 

1 V A s/ IjJ 



q -by-q matrix 

q -l" by_q 2 matrlx 


[A 1 B 

L_ 

MV • 

VV 

! (fcf 

i 


c |d 

q 2 -by-q 1 matrix 

i 

| matrix 




( 15 ) 
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a 


and P (2) 

are the usual transfer probabilities between sub- 
regions lying in the same medium, they can be denoted by 1 and 
respectively, to be consistent with the previous notation. If 
there are q^ = 6 subregions in medium I, for example, and q 9 « 4 


subregions in medium II (q = q. 


1 + *2 


= io), then equation (15) reads 







p [ 1] 

4 ir 



P l,7 

p i,e 

P l,9 

P l,10 




P ( 1} 

4 1 ’ 

P (l) 

r 3 

4 1 ' 



P 2,7 

P 2,8 

P 2,9 

P 2, 10 

1 ■ te ) 1 

pU) 

p (l) 

r 3 

p f> 

P (i) 

r 2 

P (i ) 

*0 

p(l) 

n 1 ’ 

4 1 ’ 

P (l) 

P (l) 

r l 

4 1 ’ 

pdi 

r 2 



P 3, 7 
P 4,7 

P 3,8 

P 4,8 

P 3,9 

P 4,9 

P 3, 10 
P 4, 10 


P (l) 

4 

p( i) 
3 

pU) 

2 

P U) 

p(l) 

0 

pH-) 



V 

P 5,8 

P 

5,9 

P 5,10 


4 l! 

p i 1} 

4 1 ' 

pdi 

r 2 


p(l) 

0 



_ P 6, 7 

P 6,8 

P 6,9 

P 6,10_ 


P 7,l 

P 7,2 

P 7,3 

P 7,4 

P 7,5 

P 7,6 



p J 2) 

pf ) 

4 2) 

4 2) " 

n /AN 1 

P 8,l 

P 8,2 

P 8,3 

P 8,4 

P 8,5 

P 8,6 


>-&r 

P (2) 

1 

p(2) 

0 

P (2) 

1 

P (2) 

2 

“ U) 

P 9,l 

P 9, 2 
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where the double- subscript elements 

are the 

transfer 

probabilities 

be- 

tween unequal 

subregions 

None 

of 

the elements in the B 

and 

c 


matrices 

are 

equal. 











With matrices 


and N s 

defined 

as 

before, 
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(16) 





N s (l - 

P) = N 

s 

0 

P 




(17) 


Matrix I is a q-by-q-unit matrix. Equation (l6) is used to fol- 
low neutrons from collision to collision in the configuration when 
ALj ^ ALjj. When convergence of (16) is slov, equation (17), which 

represents q simultaneous algebraic equations for the total number of 
scattering hits in each subregion, is used. 


Combined application of diffusion theory and numerical method . - 
Assume, for example, a thick slab having an isotropic source plane at the 


* 
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left extreme boundary and located in vacuum. Ordinary diffusion theory 
applies accurately in the region of about three or more mean free paths 
from the source. Denote this as region II. The entire slab is solved 
by ordinary diffusion theory. The resulting flux, however, is plotted 
only in region II because an exact transport-theory calculation would 
yield substantially the same curve only in this region. This leaves the 
region from the source plane up to three mean free paths to be worked by 
higher order theory. 

The flux in region II includes neutrons that have crossed and re- 
crossed the interface any number of times. The angular distribution 
into region I arises from this flux and is calculated by equation (B3). 

It therefore represents the final angular distribution and is consequently 
calculated only once. Thereafter, region II is treated as a vacuum in 
following the neutrons in region I. 

The first-collision distribution in region I, due to the angular 
distribution at the interface, is obtained by equation (B6). To this is 
added the contribution from the original source plane at the left bound- 
ary. By knowing the first-collision distribution, successive-collision 
densities are calculated in the usual manner, utilizing the transfer 
probabilities of region I. The resultant numerical flux completes the 
distribution for the entire slab. 

Calculation of various quantities of Interest in a slab calcula- 
tion. - From the total number of scattering hits n s (i) in each subregion 
of the given slab, the average total flux is computed by equation (6). 

The total absorption rate in each subregion and consequently in the 
entire slab is computed by 


*s 

^capt 



n s (i) 


(18) 


Since the number of neutrons per second entering the given slab - 
that is, the production - is known, then 

Total leakage = Production - Absorption 
where the absorption term is given by expression (18). 

The total number of neutrons n T transmitted outward through the 
extreme left or right boundaries from one or more collisions within the 
slab is 
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£ n s (i)pT( ti ) , ^ 


i=l 


Z = ng 


1=1 


where l r = ^r - l^AL and P T (Z) Is defined by equation (A6). 

The angular distribution of these neutrons is (appendix B)i 

| ^hl ) d P = S nS (l) e ^ 


i=l 

q. 


1 dp 

2 ^ 


P h < 0 




dp 


P h > 0 


i=l 


J 


(19) 


where | p^ | is given by equation (Bl). 

If neutrons have been followed from collision to collision according 
to equation (4), n^(i) is known for each collision. The number of neu- 
trons transmitted after exactly K collisions follows by replacing 
n s (i) in equations (19) by n|(i) and n T by n^. 

Boundary conditions at extreme bo undary of a cell. - Consider an 
infinite repetitive -slab array consisting of identical units called 
cells. At the extreme boundary of a cell, tie leakage from the cell in 
any direction is compensated by leakage into the cell in the opposite 
direction. This necessitates determining the angular distribution at 
the extreme boundaries from the flux within (eqs. (B2) and (B3)) and then 
reflecting the angular distribution. The uncollided source neutrons must 
be added to n(x-p, |p-^|)dp to complete the angular distribution. The 

first-collision distribution of these neutror s is then calculated by 
applying the technique of equation (B5) , and so forth. 


NUMERICAL EXAMPLES 
Description of Slabs Studied 

To illustrate the preceding development, the slab of sketch (a) was 
chosen with a = b = 1 centimeter. The total mean free path and the 
scattering probability per collision were respectively taken as 
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A 1 = A 11 = 0.78125 and. ( A/A s ) X = (A/Ag ) 11 = 8/l0. Each slab was divided 
into ten equal subregions (q^ = q^ = 10) giving AL^ = AL^ = 0.1280* A 
unit isotropic source plane was chosen at the interface Xq. 

In the general case where ALq ^ ALg, this conf iguration might he 
treated as a two-slab problem (each half separately) by using the angular 
distribution at xq to pass from one slab into the other* However, in 
the present case, there is the option of a single-slab treatment, which 
is, of course, much simpler* Both procedures were adopted for this 
example and gave identical results as anticipated. For the double-slab 
case a grid of 50 solid angles was used to pass from one slab into the 
other. As discussed in the text, transfer probabilities between sub- 
regions of unequal width can always be used as an alternative to the 
angular-distribution approach. 

Let P^ and P^ be defined as the respective probabilities that 
a source neutron will be transmitted through the right or left boundaries 
after exactly K collisions. For example, would be the probability 

of a transmission through the right boundary directly from the source. 
From symmetry, 



The multiple integrals for P^ and P^ are formulated in appendix C. 

T 

By employing the previously given values of A and P 0 and 

P p were evaluated as P^ = 50.0 and P^ = 47.07 per thousand source 
neutrons. Part of the contribution to Pq consists of neutrons that 
suffer the single scattering collision to the left of the source plane. 
The fraction of these that escape through the right boundary P^ was 
evaluated as 12.3 per thousand source neutrons. In addition to the 
neutron flux, these values were compared with the results of the numeri- 
cal method. 

The block diagram for the numerical analysis, using the single-slab 
treatment, is given as appendix D. Twenty collisions proved sufficient 
to give a converged flux. 

The neutron flux and transmission probabilities P^, P^, P^ T were 
obtained by a Monte Carlo calculation (appendix E) considering 10,700 
histories for up to 20 collisions each. The results were compared with 
those of the numerical method. A P3 spherical harmonics solution 
(appendix F) and a Pq diffusion-theory calculation were also performed 
for the total flux and were compared with the other methods. 
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The slab half-thickness was successively changed to l/5 and 5 total 
mean free paths with the same values of A and A/A s . The numerical- 
method flux was compared with diffusion theory for these cases. 

As a second example, the configuration cf sketch (a) was chosen, 
but with 


and with the source plane at instead of x 0 . The widths of slabs 

I and II were taken as 2.0 and 2.0177 centimeters, respectively. By 
choosing q^ = 34 subregions in slab I and q^ = 10 subregions in slab 
II, equation (ll) is satisfied. Thus, ALj = ALjj = 0.058824. There are, 
of course, many other integral pairs q^ and q^ that could have been 
chosen to achieve the equality of the Aids. 

The set of transfer probabilities Pq, I]_, P^, . . ., P44 was cal- 
culated from equation (A8a). The additional set p^, p^, . . * , p^ 

was calculated from (A8b) in order to determine the first-collision dis- 
tribution of the neutrons emanating from the source plane at X^. The 
successive-collisions technique of equation (13) was applied for 28 
collisions and yielded converged results. The number of neutrons trans- 
mitted after exactly K collisions past X^ and Xr was recorded in 
the process. 

For comparison, a Monte Carlo calculation was performed for the 
total flux and the number of transmissions. dp to 28 collisions were 
allowed, and 25,000 neutron histories were followed. 

Finally, the flux was obtained by solving the 44 simultaneous equa- 
tions represented by equation (14). 


Results 

The results of the numerical and Monte Carlo methods for the symmet- 
rical slab of sketch (a) are contained in tables I and II. Table I also 
includes results of the P^ approximation. The fluxes are plotted in 
figures 1 and 2. Figure 2 includes a P4 approximation for comparison 
with P3. 

T T 

Against the analytic values of Pq = 50.0, P^_ = 47.07, and 
P^ = 12.3 per thousand source neutrons, Monte Carlo yielded averaged 

values of 50.0, 47.81, and 11.97 per thousand as compared with 49.6, 
47.03, and 12.18 by the numerical method. Tie Monte Carlo flux and 
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transmissions sometimes suffer an appreciable deviation from symmetry in 
each half of the slab (table II), but the averaged values are in good 
agreement with those obtained by the numerical method, which, of course, 
suffers no deviations from symmetry. The total leakage was 0.43406, 
0.43556, and 0.43598 per source neutron per square centimeter per second 
by the Monte Carlo, numerical, and methods. The total absorptions, 

satisfying neutron conservation within 0.04 percent, were 0.56569, 

0.56402, and 0.56402, respectively. 

In addition to the P-]_ and P 3 fluxes, figure 2 contains the re- 
sults of the numerical -method solution utilizing 100 subregions instead 
of 20, as presented in figure 1. The finer mesh was chosen to test the 
accuracy of lumping the scattered neutrons at the center of each sub- 
region and to acquire a more detailed variation of flux than appears in 
figure 1. When averaged over each subregion of figure 1, the numerical 
flux differs by a small fraction of a percent from the 20-subdivision 
case, Implying a negligible lumping error. The same is true for the 

pT*s, K = 0,1, . . ., 19. Thus, for practical purposes the 100-subregion 
K 

case represents a precise solution of integral equation (10) and hence 
the transport equation, in the sense that further subdivisions would yield 
negligible changes and that the numerical method has exactly treated the 
angular dependence. 

Utilizing the same total and scattering mean free paths as before, 
figures 3 and 4 contain numerical -method flux plots for slabs of lj 5 and 
5 total mean free paths half-thickness, respectively, from the source • 
plane to either boundary. The diffusion- theory flux is plotted for 
comparison. 

Figure 3 shows the agreement between the numerical method and dif- 
fusion theory to be poor everywhere, whereas in figure 4 the agreement 

is excellent beyond 2ij? mean free paths from the source. 

A total of 10,700 Monte Carlo neutron histories required about 2 
hours of IBM 653 machine operating time as compared with 3 minutes for 
the numerical method. Sketch (a), worked as a double -slab problem, 
required about 4 minutes per slab per interface crossing. Twenty colli- 
sions were followed within each slab. Thirty minutes machine time corre- 
sponding to eight interface crossings yielded results identical to those 
listed in table I. 

The numerical, and Monte Carlo results for the second example (two- 
region unsymmetrical slab) are listed respectively in tables III and IV, 
and the fluxes are plotted in figure 5. Table III shows that 28 colli- 
sions proved sufficient to obtain the neutron flux accurate to the fourth 
decimal place in comparison with the solution of the 44 simultaneous 
equations for the neutron flux. 
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Figure 5 shows that the numerical flux pLots into a smooth curve 
about which the Monte Carlo solution oscillates. The reasonably close 
agreement implies a check on the correctness )f the numerical solution. 
Identical remarks apply to the transmission probabilities in tables III 
and IV. 

The numerical -method leakage through the X L and X R boundaries 
was 0.68969 and 0.061656 per source neutron to give a total-leakage 
probability of 0. 75135 per source neutron. The respective Monte Carlo 
values are 0.694384, 0.060358, and 0.75474. r ?he total absorption per 
source neutron is 0.248628 by the numerical method and 0.245178 by the 
Monte Carlo calculation and thus checks neutron conservation. 

Following 28 successive collisions numerically required about 30 
minutes machine operating time as compared with 15 minutes machine time 
in solving the 44 simultaneous equations. Twenty-five thousand Monte 
Carlo histories required about 10 hours. 


Discussion of Resull s 

The 100-subregion case of the slab of sketch (a), as opposed to the 
20-subregion case, yielded virtually identical results for all the flux 
Q*frd transmission values. From this, it is cor eluded that the error 
introduced by lumping the scattered neutrons st the center of a subregion 
was negligible for this problem. 

Anticipated qualitative features are confirmed in figures 2, 3, and 
4. The readily obtained numerical result is in excellent agreement 'with 
P 3 beyond 1 mean free path from the source and with diffusion theory 
beyond 2.5 mean free paths. From the source to about l/5 of a mean free 
path, the agreement of P 3 and P-,_ with the numerical method is poor 
and becomes worse, in general, as the source plane is approached. In 
this vicinity, a proper treatment by transport theory requires a large 
number of spherical harmonic terms to account for the large forward bias 
in the net neutron current. 

The quantity l/A = 2 was chosen such that Pq = 50. 0/l000^ and there- 
fore 2=1. 280. This value of 2 was obtained by visual interpolation 
from a set of curves given in reference 5. 

The simultaneous-equation approach of equation (14) proved best in 
the unsymmetrical-slab case or second example. This solution represents 
the flux that would be obtained by following ai infinite number of 


E-170 



E-170 


23 


successive collisions according to equation (13). This approach, there- 
fore, represents the solution to the problem. 


Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio, December 11, 1958 
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APPENDIX A 


CALCULATION OF TRANSFER PROBABUTT iES AND LUMPING ERROR 

A neutron is scattered or torn isotropically at x' in dx' in 
sketch (t). Denote by P^(x',x) the probability that it will pass 
directly through the infinite plane at x without suffering a collision. 



Sketch (t) 


The probability that it scatters with c.irection cosine p in dp 
is i dp (appendix E). The probability of traversing a distance p with- 
out a collision is e" p /\ Summing over-all angles to the right of the 
source gives P T ( x ' , x) . 

P T (x',x) = i f dp e“( x-x ' x > x' (Al) 

La */0 

The probability that the next collision occurs in dx about x is 

P(x',x)dx= lim P T (x',x) - P- r (x' x+Ax) 

Ax -*-0 

= lim | f 1 dp e _ ( x_x ' ‘/V (l - e-Ax/V) 

Ax -*• 0 2 

= 1 f 1 dp e -( x ^ x ')/V £:c 

2 Jo ^ 1 


n; t- 
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By letting t = - and. I = - ~ X , this expression becomes 

A 


P(x' ,x)dx = 


and (Al) becomes 


2A J(x-x')/A 


P T (^) = 4 


dt x > x' 


■U 


Equation (A2) represents the kernel of integral equation (10) when 
x > x' and is the well-known exponential integral. If x < x', the only 
change is to replace x - x' by x' - x. The symbol I is the number 
of mean free paths in |x - x' | centimeters. 

Integrating by parts gives 


P T h>-KM [ *r 


The value of the second term is given in Jahnke and Embde as 
1. 92X10”®. Integrating this series term by term gives 


I-L = 1.92X10" 8 + In 15 + X) ^ ' nin ' 1 " In 2 + 

n=l 


y' (-D n+1 z n 

^ n’n 


The first three terms are lumped into a single constant. The result is 
-0. 57721560. Equation (A4) becomes 


n+l,n+l 


P T (2) = 0. 5 - 0.21139220 l + -k l In l 


-u“ T n 


2 (n+l).'n 


This expression is the prohahility that a neutron scattered or horn 
isotropically at x T is transmitted through a boundary l mean free 
paths away with no intervening collisions. 
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By referring to sketch (c), the transfer probability between the 
two shaded subregions is computed as follows: 

All neutrons scattered in AL' are luaped at x* inside AL'; x* 
is ordinarily chosen at the center of AL’. This is not an essential 
choice, but it is a convenient one. The subregions have been taken suf- 
ficiently small so that the lumping error can be neglected. 

The quantity L’ is the distance in mean free paths from the source 
plane to the center of AL; L r is the distance from the center of AL’ 
to the center of AL. The probability that a neutron scattered in AL' 
will make its next collision anywhere in AL is given by 

p(iU) = p T (r) - ? T (z) (A7) 

But l ' = L' - and l = L' + Substituting from equation (A6) , 

(A7) becomes 

P(L',AL) = 0. 21139220 AL+-|^L' - y^ln^L' + ^r) ln^L' + ^ + 



p(L l j AT.) is the transfer probability from a source plane into a subregion 
of width AL whose center lies L' mean free paths away from the source 
plane. 

The source plane is ordinarily taken at the center of a subregion. 

Then, 


n; t 
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Transfer Probabilities Between Unequal Subregions 

Suppose that AL 1 = ALj lies in medium I and AL = ALjj lies in 
medium II. It can be proven that, for slab geometry, the transfer proba- 
bility still depends on the total distance in mean free paths from the 
source plane in AT, j to the collision center of ALjj*. From sketch (d), 


Slab I Interface -7 Slab II 



x 0 

Sketch (d) 
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L' = Z^ + Zp 
l T ■ ( r l ' |) fiL I 

! ? - ( r 2 - 2)^11 

Equation (AQ) now reads : 

PCAL^ALjP = 0.21139220 AL n + \ (l J + Z* 1 - ^p)ln(lj + zj 1 - ^p) - 

I ( ir i + ,11.^1)41 +l n + ^) + 





where l £ and are given by (A10). Equation (All) is used to calcu- 

late transfer probabilities between unequal subregions. 


Lumping Error 

Again, in sketch (c) the lumping error per collision e arises 
because neutrons scattered anywhere in AL 1 .are lumped at a source plane 
at the center of AL r . If S(L ! )dL f is the :rue source distribution 
inside AL ! , where dL 1 is of infinitesimal width, the lumping error for 
the transfer probability P(L r ,AL) is 

LL ! 

P(L r ,£L) (A12) 


Equation (A12) implies that if all neutrons inside AL T were actu- 
ally scattered from a source plane at its cenner, the lumping error for 
P(L r ,AL) would be zero. This is verified by substituting 
S(L T ) = S6(L ’ - L r ) into (A12). 

The maximum possible lumping error that could arise would occur if 
all neutrons were scattered in AL T at either of the two boundaries of 
AL 1 . Then, 
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S(L')dL' = SS\L' - L r - ^JdL’ (for e a ) 

or 

S(L')dL' = S5^L’ - L r + ^dL’ (for e b ) 

Then 

e a = |p(L r + 4 £,Al) - P(Lr,AL) 

€ b = |p(L r - ^,Al) - P(L r ,AL) 

For Pq, for example, the maximum possible error is 

W - P 0 - Pi- =1- (A9) - =!• (A8t) « f in 2 + g ‘ f) 

In actual practice the error will be considerably smaller than this 
because neutrons will be distributed throughout AL 1 , 

Insight into the magnitude of the lumping error may be gained by 
writing S(L t ) = £L ! + rj in equation (A12). This accurately approxi- 
mates the correct but unspecified source distribution in AL 1 , provided 
AT. * is small. The quantity r) is chosen so as to normalize the source 
distribution in AL', while £ remains an unknown parameter unless the 
source distribution is specified. The integration in (A12) can then be 
performed with the result being indicated as f(L,AL,£). Thus, 

e(£,L r ,AL) = f(L r ,AL,0 - P(Lg,AL) 

where Lg has "been written for L r in the P(L r ,AL) term of equation 
(A12). By setting e = 0 and solving this complicated equation for Lg, 
the exact location within AL' where all the scattered neutrons can he 
lumped with zero error is obtained for a given L, and AL. This 
point will rarely coincide with the center. 
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APPENDIX B 


ANGULAR DISTRIBUTION AND FIRST- COLLISION DISTRIBUTION 
DUE TO AN ANISOTROPIC SOURCE PLANE 
Angular Distribution 
Slab I Slab II 



The angular distribution at the interface of neutrons entering slab 
II from slab I is computed from the known nunber of scattering hits 
n s (i) in slab I. 

In sketch (e), is a unit vector in Ihe direction of the neutron 

velocity that is specified by the cosine of the angle between and 

the x axis; that is, x. The vectoi ^ can be considered as 

the radius vector of a unit sphere about the point P as its center* 

The surface of the hemisphere to the ri^ht of the interface is ^ 
divided into a grid of g strips of equal aieas formed by rotating 
about the x axis for h = 1,2, . . . , g. Ihen, 


i 


^0 


= 0 


Kl = K-i| + - \ ^,1) 

h = 1,2, . . 


(Bl) 


g 
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where 6h 1 i s the Kronecker delta. Thus, if g = 50, for example, 
dp = 0.02 and p^_ = 0.01, pg = 0.03, p^ = 0.05, ... = 0.99. 

The fraction, of neutrons scattered into dp is g dp. The fraction 
that reaches the interface from the scattering center in subregion i 
without further collisions and with direction cosine | p^, | is 


Thus, 


n(x Q , |p h | ) d P ^(i) J >3-P e ^ 

i=l 

>q 1+ l-i » [(4l * 1 - h - |]^I 


Similarly, 


-I 11 

n(x Q , I Phi ) dM - = nS (i) ^ d-P e <ll ^^ h 

i=q L +l 


7 11 

i-ll 


[(i - q x ) - 


^h > 0 


Ph < 0 


(B2) 


(B3) 


Equation (B2) is the angular distribution entering the second slab 
due to the flux in the first, and (B3) is the converse. It should be 
noted that the angular distribution is not isotropic. 


First-Collision Distribution in Second Slab Due to Angular Distribution 


^h 


Define 6? as the probability that a neutron with direction cosine 
will pass through a single subregion of slab I without a collision. 


-AL/ | p h | 

s 6 

e-^Il/IPhl 


(B4) 


The quantity 1 - &u is the probability that a neutron with direction 
cosine p h will collide within a single subregion. 

The number scattered from the first collision in each subregion i 
of the second medium is therefore 
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g 

n l^) = (AY 1 n(x 0 ,n h )dn(&J I ) 1 " <l1 ' 1 (1 - bj 1 ) p h > ° (B5) 

V s/ h=l 

i = q.-^+l , q-^+ 2 , . . • , q 
From slab II into slab I, 

I S 

n ! (i) = ( t ) S n ( x 0,|^hl)^(6j) <l1 ’ (1 - &£) < 0 ( E6 ) 

VS/ h=l 

i = 1,2, . . . , q. ± 

It -i T 

The factor (&^) (l - 6j^) gives the fraction of the n(xQ, inter- 

face neutrons that pass through q^ - i subregions of medium I without 
a collision and that collide in the next subregion. 
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APPENDIX C 

EVALUATION OF INTEGRALS FOR PROBABILITY OF TRANSMISSION 
WITH NO COLLISION AND ONE COLLISION 

The quantity is equivalent to equation (A3) . Fo^ the slat) of 

sketch ( a ) , l = a /A = l/A, since a has "been chosen to be 1 centimeter; 
Pq was arbitrarily chosen as 50.0/1000; A was determined so as to give 
this value, and turned out to be approximately 0.78125 (ref. 5). 



To evaluate P^, consider sketch (f). Let P^® and denote 

the respective probabilities of a source neutron suffering its first 
collision to the right or left of the origin and then escaping through 
the right-hand boundary. Then 

p^ = P*® + P?j® (Cl) 

To evaluate P*®, for example, write the expression for the number 
of neutrons per square centimeter per second (N^'(Wq,W^,P0j)Pi) ( ^0 j ’ ( ^ 1- ,< ^ (:) 0^ 
born into solid angle <T7 q, scattered into dil-^ at a radius pq from 
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the origin, and traveling the remaining distance p-^ to the right-hand 
boundary without a subsequent collision. If a total of N neutrons per 
square centimeter per second are isotropically born from the source, 

W T (fi 0 ,ft 1 ,p 0 ,p 1 )dn 0 dfl 1 dp 0 = N e“P:/ A ^ e"Pl/ A 


a - X-, 


where dft = - dp, da, p n = — L , p. =■ 

' M 0 pi l 

up all Uq t s, Hj_ r s, and x^’s, is obtained: 


, anl \i - cos | . By summing 


ff 

x-^-0 J M-q— 0 JcLq-0 
Similarly, 



dc^dpQdxi 



a-x-, 


4ttA s p 0 


_1_ 

4:rt 


A 


da^dp^ e 1 


a^o 



-i r 2* 


fl 


a-x 


1 


1 e W 
da 0 dp 0 d Xl — — — 

Xf=-a 0 JcLq=0 «/a-j = 0 


1 Pl A 
da l d Pl 4^ 


Integrating these expressions over cv, cl , ard x pives 

0 1 1 ^ 


,2 


n v 


a 


-:/v -l/v N 
e x - e u 


^ ' 4AAs /, I dv ° dVl r^r 


WC2) 


pLR = 


1 4AA F 


•A/a /*A/r 


'0 Jo 


l/\ 1 

dv 0 dvi e x 


("0 


1 - e 

1 + v o/ v l 


Integrals (C2) were evaluated on the IBM 653 by using the definit 

of double integrals and a value of A/A s equal to 0.8. The results are: 


ion 


pRR _ 34.76 __ fl pLR 12.31 

F - tuou ~ 311(1 u - iuuo~ 
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Thus, from equation (Cl), 


The values 


P T = 47.07 
1 1000 


,T _ 50. 0 p T 47.07 pLR _ 12.51 

0 " 1000" 1 1000 J 1 1000 


vere to "be checked hy the Monte Carlo and numerical methods. 



APPENDIX D 


NUMERICAL TREATMENT OF SLAB OF SKETCH (a) 

The values of the constants (see p. 12) are q = 20, a =1 centimeter, 
A - 0. 78125, A/A s = 0.8, and S x ^ = 1 neutron per square centimeter per 

second. Thus, AL = 0.1280. In evaluating the transfer probabilities, 
fifteen terms of the series in equations (A8) and (A9) were used, so 
that extreme accuracy was obtained for AL = 0. 1280. 


The accompanying block diagram summarizes the main points of the 
illustrative slab calculation but omits the steps for obtaining pLR # 

From r^(i) = n s (i), the average flux in i was calculated by 1 

equation (6). 
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NUMERICAL -METHOD BLOCK DIAGRAM 


J 



Sixty locations are used for this purpose: 20 locations to 

accumulate the number of scattering hits in each subregion for all 
collisions, 20 to store the number from the previous collision, 
n|(i), and 20 as the working area for n^ +1 (i). 
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APPENDIX E 


MONTE CARLO CALCULATION OF NEUTRON FLUX nND TRANSMISSION INTEGRALS 

Sketch (g) shows one of the symmetrical, halves of the slat of sketch 
(a). The purpose of the Monte Carlo calculation was to obtain the total 



neutron flux and transmission integrals P^, pT, and P^ (whose values 

are given by eq. (C3)) and to compare them with the values and time taken 
by the numerical method. The probabilities p£, K = 2,3, . . .,19 were 
also recorded and compared, though not analogically evaluated because of 
the complexities. The fact that these multiple-scattering transmission 
integrals are so readily obtainable by Monte Carlo demonstrates the 
utility of following individual neutron hisl ories for thin slabs. 

The complete slab was divided into 20 subregions, and the total 
number of scattering hits (proportional to the average flux) was recorded 
in each. In following individual histories, all collisions were treated 
as scatterings by the well-known technique cf weighting the colliding 
neutron by (l - ZqJTL) per collision. A neutron leaving the slab from 
the K^ 1 collision was recorded in the T k cr Tk location of the 
machine, depending on the boundary through vhich it escaped. 
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To achieve accurate results, 10,700 histories were followed. 1 By 
averaging corresponding fluxes and transmissions in both halves, the 
equivalent of 10,700 case histories in each half is obtained. A Monte 
Carlo treatment of this slab involves the following: 

Choosing random numbers : 

A pseudo random number Rn-f 1 was generated by taking Rn + l*pRn 
where p = 9,677,214,091 and R 0 = 6,250,739,481 and by retaining only 
the right ten digits in the product. The extreme right digit always 
ends with a 1, but nine-digit accuracy is more than sufficient. The 1 
prevents degeneracy from occurring too quickly. 

IhL Choosing direction of travel after isotropic scattering colli- 
sion in laboratory system: 


Because the direction of travel after an isotropic scattering colli- 
sion is unrelated to any previous direction, the x-axis can always be 
chosen as the fixed reference or initial direction. The scattering angle 
thereby formed is denoted by £ and its cosine by p. 

From sketch (h), the probability distribution function for a neutron 
passing through the shaded area dA is 

p(|)d| = dA/4it 

= ^ sin | d£ 



The cumulative distribution function is 
equal to a random number Rq to give 


fo paoas’ 


and is set 


cos | = p = 1 - 2Rq 


(El) 


pThe standard deviation of P^, for example, is 
a = V n P(! " P) = yi/10,700 rj 0. 05X0^95- Thus, in the Monte Carlo calcu- 

latior Pn = — ± — ^1^2. = 50* PUUIl with 50-percent expectancy. 

u 1000 10,700 1000 
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Thus, to pick from an isotropic distribution in the laboratory system, a 
random number is chosen, and equation (El) Is applied. Picking an 
aximuthal angle is unnecessary for this slab problem. 

ll L Picking distance of travel to next collision ; 

The probability distribution function for a neutron traveling a 
distance p without a collision and then colliding in dp is 

p(p)^P = e The cumulative distribution function is p(p T )dp', 

and equating to a random number gives 

p = -A In (1 - R') = -A In R 2 (E2) 

(4) From sketch (g), the penetration distance is ; 

^K+l = X K + Pk+1 11 K 

where X K is the X coordinate of the previous or K^ 1 collision. 

The Monte Carlo problem was programmed according to the block 
diagram given on the next page. After following 10,700 histories for 
up to 20 collisions per history, the total number of scattering hits 
n^(i) in each subregion i was punched out in addition to the trans- 
mission probabilities P^, P^, P^, K = 0,'., . . ., 19. Then, the 
average flux in i was calculated by 

$(i) = n s (i)/(z s AXx:.0,700) 

where AX is the thickness of subregion i in centimeters and 10,700 
is the normalizing factor. The number of neutrons transmitted was 
also normalized to one source neutron by dividing by 10,700. 



CY-6 E-1T0 
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MONTE CABLO BLOCK DIAGRAM 
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APPENDIX F 

P 3 SPHERICAL HARMONICS SOLUTION 

The P 3 approximation in the Legendre series expansion of the 
transport flux. 


F(x,(i) 



F l (x)P l (p) 


F Z (x) 



F(x,n)P z ( M )dp 


gives rise to four simultaneous differential equations for the total 
neutron flux Fq(x) when substituted into the transport equation (ref. 
1). For isotropic scattering in the laboratory system and a distributed 
constant isotropic source S, these equations are: 


F' (x) + h Q F 0 (x) - S = 0 


Fq(x) + 2F£(x) + Sb-jF^x) = 0 
2Fj_(x) + 3F 3 (x) + Sh^F^x) = 0 
3Fg(x) + Yb^^Cx) = 0 




J 


(FI) 


The net neutron diffusion current is Fq(x), while F 2 (x) and F 3 (x) are 
higher order fluxes and currents; bq is the macroscopic total cross 
section and bQ the macroscopic absorption cross section. An identical 
set of equations holds for another region, say A, except that 
neutrons per cubic centimeter per second replaces S and the b's are 
replaced by a's. 

Consider sketch (i). As discussed in reference 1, the source plane 
is treated as a vacuum region of width 2Xq, which is eventually 
allowed to approach zero with 2xqS^ approaching 1 neutron per square 
centimeter per second. Since region A is vacuum, a Q = aq = 0. 


E-170 



CY-6 back E-170 
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Vacuum 


(l Neutron/cc/ 


sec 


A 


a Q - 0 
a n = 0 


V 


Source 

region 

A 


-*8 x i 


I 


Region B 


b q = 0. 2560 
b x = 1.280 


Vacuum 


Sketch (i). - Slab of sketch (a) by spherical harmonics. 


The solution of equations (Fl) for regions A and B is: 

Ifc*) - ^ ' - Mx| 


F^x) = S x + A 
l v A 1 


Po(x) = A. 


F^(x) = -(2/3) S A x + A. 


A 3 


= B v e 

V=1 


F = N s y B e 

1 Z-» v v 

v=l 


F 2 = 2 8 vV 

v=l 


0 V I X 


M x ! 


F 3 = 


V=1 


. B e 
v v 


0vl x 


>(F2) 


where A and B are arbitrary constants. Also 
b. 


r = - -2 
Tv Pv 


3b A b. 

6 = 0 1 _ 1 

v 2P 2 2 

*v 


14b 1 rv 


3b 0 b l\ 


4 - "lPv + *2 » "1 - 3b 0 b l + f b l + f b 0 b i; *2 = ¥ b 0 b l 


(F3) 


The boundary conditions are (ref. 6): 
At x = 0, 
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At x = x-p 


At x = x^. 


f£(o) = 0 

F^(0) = 0 


®t( x l) = ^( Xl ) m = 0,1, 2, 3 

F ( x 2 ^t i ) p i(l a ) d 4 = 0 

f F ( Xg , ja ) P 3 ( p. ) d+o. = 0 

* / 0 


(F4) 


(F5) 


(F6) 


Equation (F4) is implied from the requirement that no net flow of 
neutrons occurs in any direction across the plane of symmetry. Equation 
(F5) expresses the continuity of the transport flux in any direction 
across the interface, and (F6) is the set of boundary conditions proposed 
hy Marshal: (in ref. l) to approximate the rigorous boundary condition 
that no neutrons reenter the slab from the vacuum; that is, 

F ( Xr> , p ) = C ),p <0. In addition, the synmetry condition F(x,p) = F(-x, -p) 
must be satisfied. 


From equation (F3), 


P l,2,3,4 



From equation (F4), 


(F7) 


A-i — 0, A^ — 0 

-L i_> 

Equation (F7) implies 


p 2 = -p x 

h = ~h 

r 2 ~ 


a 2 = 

-A 

6 2 = 

8 4 = 8 3 

U = 

" r s 

>• 

ll 

-A. 


Boundary conditions (F5) and (F6) lead t) the following set of 
equations for the arbitrary constants: 


(F8) 


fcd 


OAT- 



E-170 
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v K ! x n 

Lj B v e = A 0 

V=1 

2 r v B v ^ v,Xl1 - s A xx 

V=1 

4 


S 5 vV 


Pv I XI 


V=1 


A x ll _ 2 


v=l 


V V 


- - 3 S A*1 


$2 (4 + 56„ - 8r v )B v e C '' , l X2 ' 


= 0 


v=l 

4 


£ (1 - 5&v + 8A v) B v e 
V=1 


M x 2l 


= 0 


y (F9) 


Equations (F9) were solved for Aq, B^, B 2 , B 3 , -^4 ^ erms 

the other known quantities. In the process, e p l x i| was approximated 
by 1 + |3|x^|, since the latter is the asymptotic expression approached 
as x^ -+ 0. The limit was taken as x^ -* 0 to give 2x^3^ = 1. The 
constants A q and A 2 were eliminated, and B^, B 2 , B 3 , and B 4 were 
finally solved in terms of b Q and bj_ and then substituted back into 
(F2). The values are: 

B x = -0.00376632 B 3 = -0.0694369 

B 2 = 0.932918 B 4 = 1.399507 

The total flux Fq(x) was plotted against x (fig. 2), and the 
total number of absorptions and leakage from the slab were computed and 
compared with the values obtained by the numerical and Monte Carlo 
methods* 
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P 3 total leakage/sec = 1 - 0.56402 = 0.43598. 

Sketch (a), worked as a two-region slab, yielded identical results; using 100 sub- 
regions instead of 20 yielded results virtually indistinguishable (small fraction 
of a percent) from those listed here. 

= 12.18/1000. 
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TABLE II. - MONTE CARLO RESULTS FOR SLAB OF SKETCH (a) FOR 10,700 HISTORIES* 

(a) Fluxes 

[Origin taken at xqJJ 


Subregion 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

- Side of origin 

2.438 

1.686 

1.414 

1.201 

0.9680 

0.3850 

0.7427 

0.6070 

0.5084 

0.3968 

+ Side of origin 

2.521 

1.796 

1.419 

1.242 

1.072 

.9065 

. 745 7 

.5941 

.5299 

.4212 

Average 

2 . 480 

1.741 

.. . 

1.417 

1.222 

1.020 

.3959 

.7442 

.6005 

.5192 

.4090 


(b) Transmissions (per thousand source neutrons) 


K 

P^xlO 3 

P^ T X10 3 

Avers ge, 

per thousand sc urce neutrons 

0 

50.47 

49.53 

50. CO 

1 

46 . 95 

48.67 

t 47 .fi 

2 

35.05 

35.47 

35.26 

3 

26 . 6b 

24.74 

25. .0 

4 

19.52 

18.37 

18.2 5 

5 

12.65 

13.32 

12.29 

6 

9.58 

8.43 

9. CO 

7 

5.74 

5.27 

5.20 

8 

3.83 

4.28 

4 . C 6 

9 

2.41 

2.99 

2.70 

10 j 

1.72 

1.93 

1.22 

11 

1.12 

.939 

1 . C 3 

12 

.770 

.719 

.'45 

13 

.493 

.478 

.^86 

14 

.374 

.345 

.360 

15 

.230 

.217 

.224 

If) 

.186 

.184 

.185 

17 

.071 

.105 

.(83 

18 

.071 

.069 

.(70 

19 

.057 

.046 

.(52 


*Total leakage/sec = 0.43406/source neutron) total absoi ption/sec * 0.56569/source neutron 
10 

(as computed from 2^ AX) . 

i-i 


^P^ - 11.97/1000. 


nJ T 




oil -a L ~zo 
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TABLE III. - NUMERICAL RESULTS FOR TWO-REGION UNSYMMETRICAL SLAB 



Transmissions 

Flux 

K 

P~ T xio 4 

P^xiO 4 

Subregion 

Successive- collision 

S imult aneous - e quat ion 





solution 

solution 





Region I 

0 

5000.00 

88.607 

1 

2.1686 

2.1686 

1 

989.99 

96.626 

2 

1.5337 

1.5337 

2 

401.50 

93.172 

3 

1.3058 

1.3059 

3 

200.73 

81.131 

4 

1.1601 

1 . 1601 

4 

112.27 

66.091 

5 

1.0522 

1.0522 

5 

67.200 

51.467 

6 

.9663 

.9663 

6 

42 . 061 

38.847 

7 

.8947 

.8947 

7 

27.149 

28.684 

8 

.8334 

.8334 

8 

17.910 

20.851 

9 

.7797 

.7797 

9 

12 . 002 

14 . 989 

10 

.7319 

.7319 

10 

8.135 

10.689 

11 

.6889 

.6889 

11 

5.559 

7.578 

12 

.6498 

.6498 

12 

3.822 

5.351 

13 

.6140 

.6140 

13 

2.639 

3.767 

14 

.5811 

.5811 

14 

1.828 

2.646 

15 

.5506 

.5506 

15 

1.269 

1.856 

16 

.5223 

.5223 

16 

.8826 

1.300 

17 

.4959 

.4959 

17 

.6146 

.9100 

18 

.4 712 

.4712 

18 

.4283 

.6366 

19 

.4481 

.4481 

19 

.2986 

.4451 

20 

.4263 

.4263 

20 

.2083 

.3111 

21 

.4058 

.4058 

21 

.1454 

.2175 

22 

.3865 

.3865 

22 

.1015 

.1520 

23 

.3683 

.3683 

23 

.0709 

.1062 

24 

.3511 

.3511 

24 

.0495 

.0742 

25 

.3348 

.3348 

25 

.0346 

.0518 

26 

.3193 

.3193 

26 

.0241 

.0362 

27 

.3047 

.3047 

27 

.0162 

.0253 

28 

.2909 

.2909 

28 

.0118 

.0177 

29 

l .2778 

.2778 




30 

.2655 

.2655 




31 

.2539 

.2539 




32 

.2430 

.2430 




33 

.2330 

.2330 




34 

.2242 

.2242 





Region 

II 




35 

0.2176 

0.2176 




36 

.2089 

.2089 




37 

.1992 

.1992 




38 

.1888 

.1888 




39 

.1780 

.1781 




40 

.1668 

.1668 




41 

.1551 

.1551 




42 

.1428 

.1428 




43 

.1297 

.1297 




44 

.1150 

.1150 
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TABLE IV. - MONTE CARLO RESULTS FOR TWO-REGION 
UNSYMMETRICAL SLAB (25,000 HISTORIES) 


Transmissions 

Flux 

K 

0 

r — 1 
X 

►d 

H 

O 

Subregion 

Results 




Region I 

0 

5088.40 

85.200 

1 

2.0935 

1 

976.320 

96.160 

2 

1.4811 

2 

382.992 

92 . 080 

3 

1.3003 

3 

190.118 

78.928 

4 

1.1233 

4 

107.404 

68.326 

5 

1.0222 

5 

69.601 

46.642 

6 

.9665 

6 

44.012 

38.743 

7 

.8733 

7 

29.416 

27.141 

8 

.8096 

8 

17.702 

20.134 

9 

.7981 

9 

12 . 614 

14.028 

10 

.76 71 

10 

8.228 

10.232 

11 

.7166 

11 

5.784 

8.857 

12 

.6441 

12 

4.009 

5.018 

13 

.5930 

13 

2.429 

3.975 

14 

.5452 

14 

1.617 

2.293 

15 

.5542 

15 

1.014 

1.947 

16 

.5142 

16 

.9172 

.9563 

17 

.4909 

17 

.4454 

.7925 

18 

.4383 

18 

.4675 

.4472 

19 

.4172 

19 

.3435 

.4270 

20 

.4530 

20 

.2552 

.2319 

21 

.3996 

21 

.2036 

.4757 

22 

.3864 

22 

.0593 

.1833 

23 

.3354 

23 

.0683 

.1395 

24 

.3418 

24 

.0599 

.0955 

25 

.3189 

25 

.0427 

.0310 

26 

.3245 

26 

.0304 

.0399 

27 

.3212 

27 

.0104 

.0448 

28 

.2792 

28 

.0055 

.0116 

29 

.2876 




30 

.2717 




31 

.2591 




32 

.2505 




33 

.2366 




34 

.2260 




Region II 




35 

0.1953 




36 

.1947 




37 

.1902 




38 

.1893 




39 

.1906 




40 

.1581 




41 

.1341 




42 

.1333 




43 

.1430 




44 

.1049 
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Figure 1. - Comparison of Monte Carlo and numerical -method solutions. 
Total mean free path, 7 \, 0.78125 centimeter per collision; scattering 
probability, X/X , 0.8. 
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Method 


Numerical ( LOO 
subdivisions ) 




0 * 2 -4 .6 .8 1.0 

Coordinate, x, cm 

Figure 2. - Comparison of P x , P 3 , and 100- subregion-numerical - 

method solutions. Total mean free path, 0.7812b centimeter 
per collision; scattering probability, h/\ , 0.8. 
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Figure 3. - Numerical and diffusion-theory curves 
for symmetrical slab of l/5 mean free path half- 
thickness. Total mean free path, A, 0.78125 


centimeter per collision; scattering probability, 

A/A s , 0.8. 
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Figure 5. - Numerical-method flux and Monte Carlo flux for unsymmetrical double-slab problem. 





